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ABSTRACT 

On December 26, 2004 a major earthquake occurred off the coast of Northern Sumatra measuring 
9.0 on the Richter scale. The earthquake triggered a Tsunami that caused more than 250,000 
casualties in the region. The historical tragedy revealed a need for urgent management for disaster 
prevention and/or relief, such as preparation of hazard maps, construction of a Tsunami warning 
system and education for local people in disaster response. 

Immediately following the event, DHI Water & Environment (DHI) simulated the propagation of 
the Tsunami wave using the global setup of the flexible mesh version of MIKE 21, with a zoom in 
to the areas where the greatest effects were reported during December 26. In the first months of 
2005 DHI, in cooperation with our Japanese partner company Alpha Hydraulic Engineering 
Consultants, have been gathering information relating to both the fault parameters of the quake and 
the initial shape of the water level displacement, furthermore, as many observations as have been 
possible. 

This paper describes the use of the additional information and the advances that has been made in 
describing the Asian Tsunami wave. 
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1. INTRODUCTION 

A Tsunami wave can typically be generated as a result of sizeable earthquakes, which generates an 
instantaneous displacement of the sea-bottom, and thus a displacement of the sea surface. This 
displacement in turn generates a wave, which behaves as a shallow surface wave. 
As Tsunamis travel across the open oceans, the amplitude of the wave of only a few centimetres is 
hardly noticeable. However, as they approach shore they will gradually build up and in areas where 
the topography of the seabed is especially suited they will reach heights of 10+ meters. The speed at 
which the Tsunamis travel of up to 800 km/h (500 mph) over the open ocean, slowing to 30-100 
km/h as they reach shore, combined with the extreme shoaling is what makes them so disastrous. 

MIKE 21 rectilinear has previously been used to simulate the Atlantic Tsunami of 1755 (ref. 
Webpage) in the waters west of Portugal. This Tsunami caused devastation to the nearby cities, 
especially Lisbon and Casablanca on the Moroccan coast. The experiences from that study, which 
was carried out as part of a feasibility and design study of a power plant on the Moroccan coast, 
facilitated the quick and optimum set-up of MIKE 21. 
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Recently MIKE 21 has also been used to study the effect of the major 1960 Chilean earthquake on 
the Chacao Channel entrance, Chile. The 1960 earthquake created a tsunami that destroyed not only 
the Chilean coast, but also damaged coastal installations in Hawaii and Japan. The tsunami formed 
up into a 15 m high bore on the Pacific Ocean beaches north and south of the Chacao Channel 
entrance. 

Following the immediate response to the Asian Tsunami disaster, DHI has analysed the Asian 
Tsunami using the state of the art modelling package MIKE 21. The rectilinear MIKE 21 version, 
as well as the flexible mesh (FM) version with parallel computing method, has been applied to 
calculate the propagation of the wave. 

The simulations have been setup with a finer geomorphological resolution in selected areas of 
interest, i.e. Banda Aceh region in Indonesia and Phuket Island in Thailand. The detailed resolution 
near land and large gradients in water depth are necessary to describe the local shoaling effect of the 
Tsunami, and especially MIKE 21 FM model makes this bathymetry modelling quite easy. 

Initial surface elevation maps for the Tsunami based on the USGS fault model, DCRC fault model 
and some original fault models have been used in the simulations. Original initial maps were 
generated from some series of fault parameters that have a magnitude of 9.1 on the Richter scale 
using the Alpha geological model with the international standard theory (Manshinha and Smylie, 
1971). Calibrations of these models are implemented with some measurements, such as those 
presented by NOAA. The effect of the land displacement has also been included in the MIKE 21 
simulations using the land slide option. 

2. HYDRODYNAMIC SET-UP OF A GLOBAL MODEL IN A FLEXIBLE MESH 

We have attempted to simulate the Tsunami propagation using an operational global 2D model with 
flexible mesh (MIKE 21 FM). With this flexible mesh we are able to zoom into interest areas where 
the Tsunami have had the most devastating impacts, in addition to the Tsunami global propagation. 
The model includes 429,459 mesh elements and the time step interval setting is 2 seconds 
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Figure 2. 1 The mesh used for the MIKE 21 FM simulation, with a zoom to various regions. 



2.1 Initial Conditions 

To create initial conditions we used the geological model (the ALPHA geological model) to 
generate initial maps from a series of fault parameters as shown in Table 2.1. 

Table 2.1 Fault parameters for initial surface elevation. 



Fault parameters 


Fault 1 


Origin X (deg.) 


1.500 


Origin Y (deg.) 


94.000 


Strike angle (deg) 


335 


Fault length (km) 


900 


Fault width (km) 


150 


Depth of the upper fault (km) 


10 


Dip angle (deg.) 


8 


Slip angle (deg.) 


110 


Dislocation (m) 


11 


Total moment magnitude (Mw) 


9.1 



The model is based on Manshinha and Smylie's theory (Mansinha and Smylie, 1971). It calculates a 
displacement of the seabed that results from seismic fault movements in the earth's crust, assuming 
that the crust consists of an elastic body and fault shape is rectangular. In this model, the 
displacement depends on Poisson's ratio. 

The initial condition for MIKE 21 FM is shown in Fig. 2.2. This gives a best and quick guess for 
the initial surface elevation, which is the most simple fault model with 9.0 Mw. 
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Figure 2.2 The initial surface elevation for FM modelling. 



2.2 MIKE 21 FM Global Results 

Using the initial field from above, the model generates a Tsunami, which propagates across the 
Indian Ocean. A few snapshots are shown below. 




Figure 2.3 Wave propagation across the Indian Ocean. 



This initial model is too coarse to simulate the Tsunami because the explicit scheme requires high 
resolution to avoid damping. Further refinements are ongoing at DHI, but we will in the remainder 
of this paper concentrate on the rectilinear version of MIKE 21. 
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3. 



HYDRODYNAMIC SET-UP REGIONAL MODEL 



We have attempted to find out an optimum initial condition for the Tsunami to explain observed 
tsunami propagation situation well with MIKE 21 (rectangular grid version). 

The model is nested simulation with 9720m and 3240m as grid scale. The time step interval is 60 
sec. 
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Figure 3. 1 Regional bathymetry. 
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Figure 3.2 Local bathymetry. 
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3.1 Initial Conditions 

We have attempted to find an optimum initial condition for the Tsunami to explain the observed 
Tsunami propagation situation well. 

First of all, we tried two initial maps presented from other institutes as preliminary simulation. One 
is from the U.S. geological survey (USGS model, 9.0Mw, Fig.3.3 Left) that results from the 
analysis of seismic waveforms. The other is from the Disaster Control Research Center in Japan 
(DCRC modeB, 9.1Mw, Fig. 3.3 Right) that results from some trials of Tsunami simulation. 




Figure 3.3 Initial surface elevation ( Left: USGS model, Right: DCRC model3 ) 

The fault parameters for the DCRC-3 fault model are shown below. 

Table 3.1 Fault parameters for initial surface elevation (DCRC model3). 



Fault parameters 


Fault 1 


Fault2 


Fault3 


Origin X (deg.) 


2.500 


5.000 


10.000 


Origin Y (deg.) 


95.750 


94.000 


92.000 


Strike angle (deg) 


329 


340 


5 


Fault lenght (km) 


330 


570 


300 


Fault width (km) 


150 


150 


150 


Depth of the upper fault (km) 


7 


7 


7 


Dip angle (deg.) 


15 


15 


15 


Slip angle (deg.) 


110 


110 


110 


Dislocation (m) 


11 


11 


11 


Total moment magnitude (Mw) 


9.1 



0 USGS model was calculated with inhomogeneous fault model. 



As next step, we tried two new initial maps based on the analysis of preliminaly results in 
cooperation with Alpha Hydlaulic Engineering Consultants Co.Ltd (ALPHA model 1 and 2, 9.1Mw 
each, Fig. 3.4, table 3.2 and 3.3). Finally, we tried the modified initial map based on the analysis of 
former simulation (ALPHA model 3, 9.1Mw, Fig. 3.4, fault parameters table 3.4). 




Figure 3.4 Initial surface elevation ( Left:ALPHA modell, Middle:ALPHA model2, Right:ALPHA model3 ) 
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Table 3.2 Fault parameters for initial surface elevation (ALPHA model 1 ). 



Fault parameters 


Fault 1 


Fault2 


Fault3 


Origin X (deg.) 


2.500 


4.800 


9.800 


Origin Y (deg.) 


95.000 


93.070 


91.730 


Strike angle (deg) 


320 


345 


5 


Fault lenght (km) 


330 


570 


370 


Fault width (km) 


150 


150 


150 


Depth of the upper fault (km) 


7 


7 


7 


Dip angle (deg.) 


15 


15 


15 


Slip angle (deg.) 


110 


110 


110 


Dislocation (m) 


15 


10 


7 


Total moment magnitude (Mw) 


9.1 



Table 3.3 Fault parameters for initial surface elevation (ALPHA model 2). 



Fault parameters 


Fault 1 


Fault2 


Fault3 


Origin X (deg.) 


2.500 


4.200 


9.960 


Origin Y (deg.) 


92.250 


93.360 


92.030 


Strike angle (deg) 


312 


347 


7 


Fault lenght (km) 


270 


640 


350 


Fault width (km) 


150 


150 


150 


Depth of the upper fault (km) 


7 


7 


7 


Dip angle (deg.) 


15 


15 


15 


Slip angle (deg.) 


110 


110 


110 


Dislocation (m) 


13.5 


11.5 


7 


Total moment magnitude (Mw) 


9.1 



Table 3.4 Fault parameters for initial surface elevation (ALPHA model 3). 



Fault parameters 


Fault 1 


Fault2 


Fault3 


Origin X (deg.) 


2.500 


5.150 


10.250 


Origin Y (deg.) 


95.100 


93.500 


92.150 


Strike angle (deg) 


315 


344 


10 


Fault lenght (km) 


320 


550 


360 


Fault width (km) 


150 


150 


150 


Depth of the upper fault (km) 


7 


7 


7 


Dip angle (deg.) 


15 


15 


15 


Slip angle (deg.) 


110 


110 


110 


Dislocation (m) 


13.5 


12 


7 


Total moment magnitude (Mw) 


9.1 



3.2 Regional Results 

We validated each result with the satellite Tsunami propagation images (2 hrs, 3:15 hrs and 8:50 hrs 
after the earthquake) from the National Oceanic and Atmospheric Administration in U.S.(NOAA). 
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Figure 3.5, Comparison between NOAA SSH, NOAA model and MIKE 21 rectilinear with initial surface 
elevations from USGS (Left) and DCRC-3 (Right). 



Preliminary results (the results for USGS model and DCRC model 3 show that the south part below 
5 deg. N Lat. of USGS result is better than that of DCRC model 3, while the north part above 5 deg. 
N Lat. of DCRC model 3 seems to be better than that of USGS model. Therefore, we suppose that 
the next trials should be similar to USGS model for south part of the map, and inversely similar to 
DCRC model 3 for north part of the map. 

For the USGS model the sea level change along the satellite trace is similar to the measurement 
around the equator at 2 hours after the earthquake. However, wave amplitude is quite lower than the 
measurement especially for northern part. That is because this initial elevation area concentrates at 
off the west coast of Northern Sumatra, and this also results in the lower wave height around the 
Phuket area. 

For the DCRC model 3 the result at 2 hours after the earthquake has a little similar trend to the 
measurement above 5 deg.N latitude. However, the wave phase is quite different between 
simulation and meaurement below 5 deg. N latitude. This difference suggests that a southern part of 
initial elevation should be located futher off from the coast of Northern Sumatra. 

The series of tests with the Alpha geological model (Alpha model 1 - model 2 - model 3) are 
calibrated to the SSH NOAA observations and model, and show considerable improvements. 
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For direct comparison to the sattelite SSH images, we have extracted the sattelite track from the 
MIKE 21 output, and compared them to the NOAA tracks. 

In the ALPHA model 1, the result describes general phase characteristics of the sea level change 
along the satellite trace well. However, wave amplitude is bigger than the measurement between 5 
deg.S and the equator. This suggests that a southern part of initial elevation should have lower wave 
height. 

Alpha model-1 Alpha model-2 

2:00 hours after Tsunami 2:00 hours after Tsunami 




Figure 3.7, Alpha model-1 (left) and model-2(right) compared to NOAA SST tracks. 
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With the initial field from ALPHA model 2, the result at 2 hours after the earthquake describes 
phase characteristics of the sea level change along the satellite trace well . However, wave 
amplitude is bigger than the measurement between 5 deg.S and the equator. This suggests that a 
southern part of initial elevation should have lower wave height. 
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Figure 3.8, Alpha model-3 compared to NOAA SST tracks. 



Finally with ALPHA model 3 initial waterlevel, the result is almost similar to ALPHA model 2, but 
this includes a little improvement of wave reproduction beneath the equator. 
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3.3 Local Results 

As described in the introduction the evaluation of the models in near shore areas requires both high 
resolution and accurate bathymetry data. We did however investigate some of the eyewitness 
accounts relayed on internet and in the media. One of these was described in EOS (see References), 
and we checked Tsunami arrival time in Sri Lanka for this result (Fig. 3.9). 




Figure 3.9 Calculated wave fluctuations around south Sri Lanka (ALPHA model3). 
(Red arrows show the Tsunami peak timing from witness in Ahungalle, 30 km north ofGalle) 



According to the witness (EOS), the Tsunami peaks were at approximately 0930, 1010(the big one), 
1110, 1 150, 1235, 1255, etc., local time at Ahungalla, Sri Lanka (about 30 km north of Galle), and 
3rd and fourth surges were larger than the first. The simulated results show good agreement with 
this witness. 



In videos recorded by various witnesses to the event one can on locations east of the fault line see 
the first arrival of the negative wave, and also see the extreme shoaling causing a continuous wave 
in exceedence of 10 meters in height, also eyewitnesses reports several waves on top of a larger 
wave. All these effects will not be captured by the hydrodynamic models, because of the extremely 
high resolution required. 

To illustrate how the height of a tsunami wave changes when it approaches the shore and inundates 
an urban area, a number of Boussinesq type model simulations have been made with DHI's MIKE 
21 BW model. This model includes the processes that determine how the tsunami (and the local 
inundation) behaves in coastal areas, for example processes like non-linear shoaling, diffraction, 
reflection, wave breaking and run-up. 
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Figure 3. 10 Wave propagation, transformation and inundation simulated by using DHI's Boussinesq Wave 
Model, MIKE 21 BW. The figures show the effect of shoaling, wave breaking and reflection from 
land and buildings. 

4. DISCUSSION 

Two different MIKE 21 model set-ups with different grids have been used to describe the spreading 
of the Asian Tsunami throughout the Indian Ocean. The propagation depends primarily on the 
method used to determine the initial conditions for the wave, and secondly on the bathymetry of the 
area. From the relatively few measurements available all indications show that both grids correctly 
describe the timing of the wave both in the deep waters and on the shallow shelf. 

The rectilinear MIKE 21 matches the amplitude of the wave in the deep ocean as compared to 
satellite determined SSH. The flexible mesh model has not been run with sufficient horizontal 
resolution to compensate for the damping that occurs with the fast travelling wave in combination 
with the explicit solution and the Courant criteria of V2, and therefore the amplitude dissipate too 
quickly on the deep ocean. 

As the wave approaches shore and the detailed local bathymetry determines the degree of 
diffraction and shoaling, which in turn determines how far inshore the wave will travel, the 
resolution of the models becomes the determining factor. Near Galle on the Sri Lankan west coast 
where the wave has been diffracted around the southern tip and the main part of the shoaling occurs 
on the west side, it appears in figure 3.9 that the model resolution is sufficient to describe these 
processes. 

For the local effects on the leeward side of Phuket in Thailand near Krabi and south thereof, the 
area consists of many islands and with a bathymetry full of details, and the resulting wave is a 
combination of the diffracting wave around Phuket with the wave that travelled down and up 
through the Malacca strait. In this case, resolution is most likely too coarse - although the fine 
resolution in the flexible mesh model seems to capture some details, and also the rectilinear model 
gets the overall picture. Unfortunately, at the time of writing we do not have sufficiently good 
measurements available to confirm the local effects in the South of Thailand. 

It has been shown that despite very limited data, the nature of the Tsunami is such that the overall 
propagation can be accurately described with the numerical engines used, as long as they have 
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sufficient resolution in space. This becomes more acute as the wave approaches shore, where the 
local bathymetry needs to be both accurately described and resolved in the models. 

It has been suggested to use models in a future warning system in this region - and the modelling 
described above confirms the feasibility of this approach. The nature of the initial elevation 
however needs, to be determined so accurately that you can not do the modelling in a fast response 
situation to anything other than the timing of impact on various shores. Instead modelling should be 
used to create an atlas of the coast describing potential high impact areas - Hazard Maps. This can 
be done by using storm surges and other events to calibrate the bathymetry of high resolution local 
models. 
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